This file analyzes the changes in S/H cell ratio (bleaching and symboint shuffling) in Pocillopora spp.: * 29 colonies were sampled in 2014 (Before ENSO), 2015 (during a 1st peak of heat) and 2016 (after a 2nd peak of heat) * These samples were collected at Uva Island, Gulf of Chiriquí, Panama.

# Libraries    
    library(plyr)
    library(tidyverse) # install.packages('tidyverse')
    # library(devtools)  # install.packages("devtools")
    # devtools::install_github("jrcunning/steponeR")
    library(steponeR)
    library(scales)
    library(gridExtra)
    library(lme4)
    library(lmerTest)
    library(multcomp)
    library(emmeans)
   
# Default ggplot settings
    ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          axis.line = element_line(colour = "black"),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

1. Symbiont to host (S/H) cell ratio and Durusdinium proportion

Get the raw data for all Pocillopora samples applying R.Cunning steponeR function

Pdam.plates <- list.files(path="data", pattern=".csv", full.names=T)
Pdam.plates
##  [1] "data/Chi-Pdam_Plate2.csv"   "data/Chi-Pdam_Plate3.csv"  
##  [3] "data/Chi-Pdam_Plate4.csv"   "data/Chi-Pdam_Plate5.csv"  
##  [5] "data/Chi-Pdam_Plate6.csv"   "data/Gal-Pdam_Plate13.csv" 
##  [7] "data/Gal-Pdam_Plate14.csv"  "data/Gal-Pdam_Plate15.csv" 
##  [9] "data/Gal-Pdam_Plate16.csv"  "data/Gal1_Dprobes_data.csv"
## [11] "data/Pan-Pdam_Plate10.csv"  "data/Pan-Pdam_Plate7.csv"  
## [13] "data/Pan-Pdam_Plate8.csv"   "data/Pan-Pdam_Plate9.csv"  
## [15] "data/Pdam_Plate11-Deep.csv"
Pdam.Out <- steponeR(files=Pdam.plates, target.ratios=c("C.Pdam", "D.Pdam"), 
                     fluor.norm=list(C=3.798, D=0, Pdam=5.416),
                     copy.number=list(C=9, D=1, Pdam=3),
                     ploidy=list(C=1, D=1, Pdam=2),
                     extract=list(C=0.813, D=0.813, Pdam=0.982))

# Target ratio results
Pdam<-Pdam.Out$result

Get colony and year information from sample labels

# Parse sample names, times, treatments, colonies, plates
    sample.year <- rbind.fill(lapply(strsplit(as.character(Pdam$Sample.Name), split="_"), 
     function(X) data.frame(t(X))))
        colnames(sample.year) <- c("Spp", "Colony","Year")
        Pdam <- cbind(sample.year, Pdam[,-1])

# Keep unique sample ID to check reruns  
    Pdam$Sample<-paste(Pdam$Year,Pdam$Colony, sep="_")
    Pdam$ID<-paste(Pdam$Sample,Pdam$File.Name, sep=".")

Calculate total S/H ratio and D/C ratio and Proportion of Durusdinium in each sample

# 0. If Clade only detected in one technical replicate, set its ratio to NA (becomes zero)
  Pdam$C.Pdam[which(Pdam$C.reps==1)] <- NA
  Pdam$D.Pdam[which(Pdam$D.reps==1)] <- NA

# 1. Rename cols and make NA=0
  colnames(Pdam)[which(colnames(Pdam) %in% c("C.Pdam", "D.Pdam" ))] <- c("C.SH", "D.SH")  

  Pdam$C.SH[is.na(Pdam$C.SH)] <- 0
  Pdam$D.SH[is.na(Pdam$D.SH)] <- 0
  
  Pdam <- Pdam[, c("Colony", "Year", "C.SH", "D.SH","Sample", "ID")]
  
# 2. Get the ratios and log 10
  # Total ratio
  Pdam$tot.SH <- Pdam$C.SH + Pdam$D.SH  
  Pdam$logTot.SH <- log10(Pdam$tot.SH ) 
  
  # C ratio
  Pdam$logC.SH <- log10(Pdam$C.SH)
    
  # D ratio
  Pdam$logD.SH <- log10(Pdam$D.SH)  
    
  Pdam$logTot.SH[which(Pdam$tot.SH==0)] <- NA
  Pdam$logC.SH[which(Pdam$C.SH==0)] <- NA
  Pdam$logD.SH[which(Pdam$D.SH==0)] <- NA
    
# Clade Proportion
  # D Proportion
  Pdam$D.Prp<-(Pdam$D.SH/Pdam$tot.SH)
  # C Proportion
  #Pdam$C.Prp<-(Pdam$C.SH/Pdam$tot.SH)
  
  hist(Pdam$D.Prp)

  #Dproportion<-Pdam[(Pdam$D.Prp<0.99 & Pdam$D.Prp>0.01),]
  #hist(Dproportion$D.Prp)

Note: Most of the colonies are dominated by one symbiont genus. Even proportions of multiple symbionts are rare.

2. qPCR Data Cleaning and quality control

 ## 1. Check and remove NTC wells
  ntc <- Pdam[which(Pdam$Colony=="NTC"), ]
  Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(ntc), ])
  
  ## 2. Check and remove NoHSratio samples
  NoHSratio <- Pdam[which(Pdam$tot.SH==0), ]
  Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(NoHSratio), ])

  ## 3. Chose bw samples ran more than once
  ReRunA <- Pdam[duplicated(Pdam$Sample),]
  
  n_RunA <- data.frame(table(Pdam$Sample))
  colnames(n_RunA)<-c("Sample","RanA")
  Pdam<-join(Pdam, n_RunA, type = "left")
  
  DuplicatesA <- Pdam[(Pdam$RanA>1),]
  
  # Remove bad replicates
  ToRem1<-read.csv("ToRemove10-28-16.csv")  # 11/1/16
  Pdam<-Pdam[!(Pdam$ID %in% ToRem1$ID),]
    
  n_RunB <- data.frame(table(Pdam$Sample))
  ReRunB <- Pdam[duplicated(Pdam$Sample),]
  colnames(n_RunB)<-c("Sample","RanB")
  Pdam<-join(Pdam, n_RunB, type = "left")
  
  # List of dupplicated samples, should have 0 rows now
  DuplicatesB <- Pdam[(Pdam$RanB>1),]
  
  # 4. Check and remove samples with high SD or late amplification
  StDe1.5 <- Pdam[which(Pdam$Pdam.CT.sd>1.5), ]
  Pdam<-Pdam[!(Pdam$ID %in% StDe1.5$ID),]
# End of qPCR data cleaning

3. Field/sample information

Get collection information and select only samples from Uva Island that were sampled in 2014, 2015 and 2016

# 1. Labeling and Factors

  ## Import treatment info from Field file and mortality / ITS2 / ORF data
  Info<-read.csv("SampleInfo.csv", header=TRUE)
  Mortality<-read.csv("MortalityData.csv", header=TRUE)
  Mortality$ITS_2<-factor(Mortality$ITS_2, level=c("C1d", "C1b-c", "D1"))


  ## Merge SH ratios and sample info
      Pdam<-join(Pdam, Info, by = "Sample", 
                type = "left", match = "all")
      #summary(Pdam)

# 2. Uva Data
      # Write how many years a colony was sampled
      Years <- data.frame(table(Pdam$Colony))
      colnames(Years)<-c("Colony","Times")
      Pdam<-join(Pdam, Years, type = "left")
  # Select data only from Uva Island    
      data.Uva<-Pdam[which(Pdam$Location=="Uva Island"), ]
  # Select colonies that were sampled three times (2014, 2015 and 2016)
      data.Uva<-data.Uva[which(data.Uva$Times>2),]
  # Remove "Ross" samples (colony in deep location (~50feet), not Uva reef) 
  # and 259 as well as 264 that were mislabeled in at least one of the sampling points   
      data.Uva<-filter(data.Uva, Colony!="Ross")
      data.Uva<-filter(data.Uva, Colony!="259-1")
      data.Uva<-filter(data.Uva, Colony!="264-1")

      data.Uva<-data.Uva[,c("Colony","Year.2","C.SH","D.SH","tot.SH",
                            "logTot.SH","D.Prp", "logC.SH", "logD.SH",
                            "Year","Sample", "Reef", "Spp", "Depht")]
      data.Uva$Year<-factor(data.Uva$Year, levels=c("14", "15", "16"))
      data.Uva$Colony<-factor(data.Uva$Colony)
      
    #summary(data.Uva)
    
  #3. Merge SH ratios and moratlity info
      data.Uva<-join(data.Uva, Mortality, by = "Sample", 
                type = "left", match = "all")
      summary(data.Uva)
##      Colony       Year.2          C.SH               D.SH        
##  256-1  : 3   Min.   :2014   Min.   :0.000000   Min.   :0.00000  
##  257-1  : 3   1st Qu.:2014   1st Qu.:0.000000   1st Qu.:0.00000  
##  260-1  : 3   Median :2015   Median :0.003005   Median :0.01442  
##  262-1  : 3   Mean   :2015   Mean   :0.019174   Mean   :0.02071  
##  278-1  : 3   3rd Qu.:2016   3rd Qu.:0.011841   3rd Qu.:0.03001  
##  284-1  : 3   Max.   :2016   Max.   :0.690021   Max.   :0.18875  
##  (Other):69                                                      
##      tot.SH          logTot.SH           D.Prp           logC.SH       
##  Min.   :0.00146   Min.   :-2.8356   Min.   :0.0000   Min.   :-5.0627  
##  1st Qu.:0.01294   1st Qu.:-1.8888   1st Qu.:0.0000   1st Qu.:-2.5001  
##  Median :0.02266   Median :-1.6448   Median :0.8620   Median :-2.0971  
##  Mean   :0.03988   Mean   :-1.6539   Mean   :0.5411   Mean   :-2.1808  
##  3rd Qu.:0.03579   3rd Qu.:-1.4463   3rd Qu.:1.0000   3rd Qu.:-1.7015  
##  Max.   :0.69002   Max.   :-0.1611   Max.   :1.0000   Max.   :-0.1611  
##                                                       NA's   :32       
##     logD.SH        Year       Sample                              Reef   
##  Min.   :-5.7262   14:29   Length:87          Coral Collection Point:36  
##  1st Qu.:-1.8381   15:29   Class :character   Shallow Millepora     :18  
##  Median :-1.6245   16:29   Mode  :character   M.intrincata Shallow  :12  
##  Mean   :-2.0059                              Ridley Rock Shallow   : 9  
##  3rd Qu.:-1.4640                              Gardineroseris city   : 6  
##  Max.   :-0.7241                              Bw Gard City and 4x5  : 3  
##  NA's   :29                                   (Other)               : 3  
##     Spp         Depht              Mortality    ID_Glynn   ID_Palacio
##  P.dam:51   Min.   : 9.00   Almost total: 7   P. cap:12   P. dam:24  
##  P.eff: 0   1st Qu.:15.00   No          :71   P. dam:27   P. eyd: 6  
##  P.ele:36   Median :17.00   Partial     : 9   P. ver:48   P. ver:57  
##  P.eyd: 0   Mean   :17.79                                            
##  P.mea: 0   3rd Qu.:20.00                                            
##  P.ver: 0   Max.   :32.00                                            
##                                                                      
##      ORF       ITS_2   
##  Type_1:69   C1d  :18  
##  Type_3:18   C1b-c:20  
##              D1   :49  
##                        
##                        
##                        
## 
      data.Uva$Mortality<-factor(data.Uva$Mortality, levels =
                              c("No", "Partial", "Almost total"))

ORF types versus morphological ID

  # Summary ID  *Pocillopora* by Glynn
  Summary_ID_Glynn<-data.Uva %>%
        group_by(Year.2, ORF, ID_Glynn) %>%
        dplyr::summarise(total.count=n()) 
  Summary_ID_Glynn
  # Summary ID  *Pocillopora* by Palacio
  Summary_ID_Glynn<-data.Uva %>%
        group_by(Year.2, ORF, ID_Palacio) %>%
        dplyr::summarise(total.count=n()) 
  Summary_ID_Glynn  

4. Community classification by genera presence/dominance

Symbiodiniaceae genera detected in the samples (Only-C, Only-D, Mixed communities)

# Add a column for the genera detected in each sample using qPCR
  data.Uva$Community<-NULL
  data.Uva$Community[which(data.Uva$D.Prp>=0)] <- "Mixed"
  data.Uva$Community[which(data.Uva$D.Prp==0)] <- "Cladocopium-only"
  data.Uva$Community[which(data.Uva$D.Prp==1)] <- "Durusdinium-only"

  data.Uva$Community<-factor(as.character(data.Uva$Community), 
                         levels=c("Cladocopium-only", "Mixed", "Durusdinium-only"))
  
# Summary genera *Pocillopora*
  Summary_Community<-data.Uva %>%
        group_by(Year.2, Community) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community
  # Summary genera detected by ORF type
  Summary_Community<-data.Uva %>%
        group_by(ORF, Year.2, Community) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community

Dominant Symbiodiniaceae genus in each sample (C-dominated, D-dominated)

In general, Pocillopora colonies are dominated by one Symbiodiniaceae type that accounts for more than 80-90 of the symbiont community. However, sometimes intermediate proportions are common when the symbiont community is transitioning from being dominated by one algal type to other (e.g. during heat stress).

For this data set, 2 of the colonies sampled in 2014 had relativelly even proportions of Cladocopium_C1bc and Durusdinium glynnii (Colony 267_2014, D proportion = 0.47; Colony 269_2014, D proportion = 0.50). This made difficult their classification into one domminant type for statistical testing and further plotting. Since these 2 colonies were later dominated by Durusdinim glynnii in 2015, we decided to include them in the D-dominated classification, eventhoug in 2014 there were not such a dominance (See figure 2a).

 # Add a column for the dominant genus
  data.Uva$DominantCommunity<-NULL
  data.Uva$DominantCommunity[which(data.Uva$D.Prp>=0.45)] <- "Durusdinium-dominated"  
  data.Uva$DominantCommunity[which(data.Uva$D.Prp<0.45)] <- "Cladocopium-dominated"
  
  data.Uva$DominantCommunity<-factor(as.character(data.Uva$DominantCommunity),
                              levels=c("Cladocopium-dominated","Durusdinium-dominated"))

  # Summary dominant genera *Pocillopora*
  Summary_Community2<-data.Uva %>%
        group_by(Year.2, DominantCommunity) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community2
  # Summary dominant genera by ORF type
  Summary_Community3<-data.Uva %>%
        group_by(ORF, Year.2, DominantCommunity) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community3

Set initial (2014) and final (2016) communities for each colony

  data.Uva$ITS_2<-factor(data.Uva$ITS_2, level=c("C1b-c", "C1d", "D1"))
  # Subset Aug 2014,  Aug 2015 and April 2016 samples
  D2014.data<-subset(data.Uva, Year.2==2014)
  D2014.data<-D2014.data[,c("Colony","D.Prp","Community","DominantCommunity", "ITS_2")]
  colnames(D2014.data)<-(c("Colony","D14.D.Prp","D14Community", "D14Dominant", "ITS2_14"))
  data.Uva<-left_join(data.Uva, D2014.data, by=c("Colony")) # Add initial community to all samples
   
  D2016.data<-subset(data.Uva, Year.2==2016) # D1a under control temperature
  D2016.data<-D2016.data[, c("Colony","D.Prp","Community","DominantCommunity", "ITS_2")]
  colnames(D2016.data)<-(c("Colony","D16.D.Prp","D16Community", "D16Dominant",  "ITS2_16"))
  data.Uva<-left_join(data.Uva, D2016.data, by=c("Colony")) # Add final community to all samples

Initial qPCR BASED analysis (No ITS2 or ORF information included)… skip to 5b?

5a. Changes in the number of Pocillopora colonies dominated by each genus

To test the statistical significance of changes in the number of colonies dominated by Durusdinium versus Cladocopium we considered inclussion and exclussion of these 2 initially even colonies, and to place them in the D-dominated or C-dominated categories in 2014.

  • Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the D-dominated category. P-value= 0.1064
Input =("
Year    C     D
2014    15   14
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.8336155 9.7658946
## sample estimates:
## odds ratio 
##   2.760908
  • Fisher test including initially (Aug 2014) even colonies. One colony was placed in the C-dominated category (D proportion=0.47), while the second one was placed in the D-dominated category (D proportion = 0.50). P-value = 0.06099
Input =("
Year    C     D
2014    16   13
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.06099
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.9556883 11.2396870
## sample estimates:
## odds ratio 
##   3.162767
  • Fisher test removing initially (Aug 2014) even colonies. N=27 colonies instead 29. P-value = 0.09778
Input =("
Year    C     D
2014    15   12
2016    8    19
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.09778
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.848568 10.659493
## sample estimates:
## odds ratio 
##   2.906792
  • Fisher test including initially (Aug 2014) even colonies. Both colonies placed in the C-dominated category (D proportion=0.47, colony 267; and 0.50, colony 269). p-value = 0.03295
Input =("
Year    C     D
2014    17   12
2016    8    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.03295
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.094473 13.006776
## sample estimates:
## odds ratio 
##   3.629777

6a. Changes in colony proportions based on detected genera (qPCR)

  • Based on detected Symbiodiniacea genera in 2014 (qPCR)
Histo_Dprop <- ggplot(data.Uva, aes (D.Prp , fill=(D14Community))) + ggthe_bw +
  geom_histogram(binwidth = 0.05, boundary=0) +
  facet_grid(~Year.2, labeller = labeller(Year.2=c(
                `2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
  scale_fill_manual(values=c("blue3","purple", "red3"),
                    guide = guide_legend(title.position = "top", nrow = 1),
                    name = "Initial symbiont community (August 2014)", 
                    labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
    scale_x_continuous(name= expression(
    italic(Durusdinium)~proportion~(D/H)/(S/H)),
                     breaks=(seq(0,1,by=0.1))) +
    scale_y_continuous((name= "Number of colonies"),
                     breaks=seq(0,15,by=5)) +
    theme(legend.position="bottom",
        legend.title.align=0.5,
        strip.background =element_rect(fill=NA))
Histo_Dprop

Traject_Proportion_Facet <- ggplot (data.Uva) +
   ggthe_bw + 
    geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+ 
    facet_wrap(~D14Community, labeller = labeller(
    D14Community=c(`Cladocopium-only` = "Cladocopium-only (2014)", 
                   `Mixed` = "Mixed (2014)", 
                   `Durusdinium-only` = "Durusdinium-only (2014)"))) +
    geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3, 
              height = 0.05, width = 0.15) +
    #labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
    theme(legend.position="none")
Traject_Proportion_Facet

Figure Not used: qPCR and colony proportion by dominant community

Histo2014_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(D14Community))) +
      ggthe_bw + coord_flip()+ ggtitle("a.")+
      geom_histogram(binwidth = 0.05, boundary=0) +
      scale_fill_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial symbiont community (August 2014)", 
                        labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
      scale_x_continuous(name= expression(
                        italic(Durusdinium)~proportion~((D/H)/(S/H))),
                        breaks=(seq(0,1,by=0.1)),
                        expand=c(0.01,0.01)) +
      scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                        limits = c(0,17),
                        breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
#Histo2014_Dprop 

Histo2016_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(D14Community))) +
      ggthe_bw + coord_flip()+ ggtitle("c.")+
      geom_histogram(binwidth = 0.05, boundary=0) +
      scale_fill_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial symbiont community (August 2014)", 
                        labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
      scale_x_continuous(name="",
                        breaks=(seq(0,1,by=0.1)),
                        expand=c(0.01,0.01)) +    
      scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                        limits = c(0,17),
                        breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
#Histo2016_Dprop

Traject_Proportion<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(D14Community), group=(Colony))) + ggthe_bw + ggtitle("b.") +
  geom_line(linetype=2) +
  geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
  
  scale_colour_manual(values=c("blue3","purple", "red3"),
                    guide = guide_legend(title.position = "top", nrow = 1),
                    name = "Initial symbiont community (August 2014)", 
                    labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
  scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
  scale_y_continuous((name= ""),
                     breaks=seq(0,1,by=0.1),
                     expand=c(0.01,0.01)) +
  
  theme(legend.position="bottom",
        legend.title.align=0.5,
        strip.background =element_rect(fill=NA))
#Traject_Proportion

Figure_2<-grid.arrange(Histo2014_Dprop,Traject_Proportion, Histo2016_Dprop,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_2qPCR.svg", plot=Figure_2, width=7, height=4)

7a. Changes in total Symbiodiniaceae abundance (bleaching): qPCR

Based on initial (2014) community

SH_2014Co<-lmerTest::lmer(logTot.SH ~ Year * D14Community + (1|Colony), data=data.Uva)
    summary(SH_2014Co)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 83.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6538 -0.5202 -0.1347  0.3700  2.9430 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.03036  0.1742  
##  Residual             0.10764  0.3281  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                     Estimate Std. Error      df t value
## (Intercept)                          -1.6402     0.1238 71.1153 -13.246
## Year15                               -0.4132     0.1547 52.0000  -2.672
## Year16                               -0.5536     0.1547 52.0000  -3.579
## D14CommunityMixed                    -0.0137     0.1638 71.1153  -0.084
## D14CommunityDurusdinium-only          0.1232     0.1805 71.1153   0.682
## Year15:D14CommunityMixed              0.7731     0.2046 52.0000   3.778
## Year16:D14CommunityMixed              0.6960     0.2046 52.0000   3.402
## Year15:D14CommunityDurusdinium-only   0.4584     0.2255 52.0000   2.033
## Year16:D14CommunityDurusdinium-only   0.3869     0.2255 52.0000   1.716
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## Year15                              0.010053 *  
## Year16                              0.000756 ***
## D14CommunityMixed                   0.933599    
## D14CommunityDurusdinium-only        0.497283    
## Year15:D14CommunityMixed            0.000408 ***
## Year16:D14CommunityMixed            0.001294 ** 
## Year15:D14CommunityDurusdinium-only 0.047152 *  
## Year16:D14CommunityDurusdinium-only 0.092125 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 D14CmM D14CD- Y15:D14CM Y16:D14CM Y15:D14CD
## Year15      -0.624                                                          
## Year16      -0.624  0.500                                                   
## D14CmmntyMx -0.756  0.472  0.472                                            
## D14CmmntyD- -0.686  0.428  0.428  0.519                                     
## Yr15:D14CmM  0.472 -0.756 -0.378 -0.624 -0.324                              
## Yr16:D14CmM  0.472 -0.378 -0.756 -0.624 -0.324  0.500                       
## Yr15:D14CD-  0.428 -0.686 -0.343 -0.324 -0.624  0.519     0.259             
## Yr16:D14CD-  0.428 -0.343 -0.686 -0.324 -0.624  0.259     0.519     0.500
    #step(SH_2014Co)
    
    SH.emmc<-emmeans(SH_2014Co, ~Year*D14Community)
    SH_groups<-multcomp::cld(SH.emmc)
    SH_groups<-SH_groups[order(
                SH_groups$Year, SH_groups$D14Community),]
    SH_groups
    plot(emmeans(SH_2014Co, ~Year|~D14Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~D14Community)

SH_2014Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(D14Community))) +  ggthe_bw +
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
      
      scale_colour_manual(values=c("blue","purple", "red"),
                            guide = guide_legend(title.position = "top", nrow = 1),
                            name = "Symbionts detectected in 2014",
                            labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.5,0.1),
            strip.background =element_rect(fill="white")) +
      scale_y_continuous(limits = c(-3, -0.8),
                name="\n Total symbiont to host cell ratio (S/H)", 
                labels = (math_format(10^.x)),
                expand = c(0.0, 0)) + 
      scale_x_discrete(name="Year")
SH_2014Com

Based on current community

SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
    summary(SH_Community)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 91.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3973 -0.5171 -0.0987  0.3905  4.0587 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.02404  0.1550  
##  Residual             0.12526  0.3539  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                      -1.60009    0.12804 76.71004 -12.497  < 2e-16
## Year15                           -0.20778    0.15540 49.73821  -1.337 0.187280
## Year16                           -0.71036    0.17979 47.25498  -3.951 0.000258
## CommunityMixed                   -0.09453    0.16829 77.96039  -0.562 0.575910
## CommunityDurusdinium-only         0.09882    0.18645 77.17622   0.530 0.597641
## Year15:CommunityMixed             0.58170    0.24003 54.72489   2.423 0.018708
## Year16:CommunityMixed             0.81742    0.24650 54.96060   3.316 0.001622
## Year15:CommunityDurusdinium-only  0.27313    0.22958 47.93526   1.190 0.240035
## Year16:CommunityDurusdinium-only  0.60594    0.24059 49.84010   2.519 0.015039
##                                     
## (Intercept)                      ***
## Year15                              
## Year16                           ***
## CommunityMixed                      
## CommunityDurusdinium-only           
## Year15:CommunityMixed            *  
## Year16:CommunityMixed            ** 
## Year15:CommunityDurusdinium-only    
## Year16:CommunityDurusdinium-only *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmntM CmmnD- Y15:CM Y16:CM Y15:CD
## Year15      -0.733                                                 
## Year16      -0.602  0.497                                          
## CommuntyMxd -0.754  0.576  0.457                                   
## CmmntyDrsd- -0.686  0.504  0.414  0.524                            
## Yr15:CmmntM  0.472 -0.658 -0.322 -0.628 -0.326                     
## Yr16:CmmntM  0.459 -0.365 -0.740 -0.608 -0.318  0.412              
## Yr15:CmmnD-  0.496 -0.677 -0.337 -0.383 -0.717  0.440  0.243       
## Yr16:CmmnD-  0.451 -0.369 -0.747 -0.326 -0.682  0.240  0.543  0.541
    #step(SH_Community)
    
    SH.emmc<-emmeans(SH_Community, ~Year*Community)
    SH_groups<-multcomp::cld(SH.emmc) 
    SH_groups<-SH_groups[order(
                SH_groups$Year, SH_groups$Community),]
    SH_groups
    #write.csv(SH_groups,"SH_groupsB.csv")
    
plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)

SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) +  ggthe_bw +
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1, 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                       position=position_dodge(width=0.95))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
      
      scale_colour_manual(values=c("blue","purple", "red"),
                            guide = guide_legend(title.position = "top", nrow = 1),
                            name = "Symbionts detectec in 2014",
                            labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.5,0.1),
            strip.background =element_rect(fill="white")) +
      scale_y_continuous(limits = c(-3, -0.8),
                name="\n Total symbiont to host cell ratio (S/H)", 
                labels = (math_format(10^.x)),
                expand = c(0.0, 0)) + 
      scale_x_discrete(name="Year")
SH_Com

anova(SH_2014Co, SH_Community)

Considering initial and current symbiont community: Figure 3: qPCR

SH_Communities<-lmerTest::lmer(logTot.SH ~ Year * D14Community * Community + (1|Colony), data=data.Uva)
    summary(SH_Communities)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community * Community + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 81.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.73407 -0.48463 -0.06579  0.35339  2.91935 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.02179  0.1476  
##  Residual             0.11026  0.3320  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                     Estimate Std. Error      df t value
## (Intercept)                          -1.6402     0.1211 70.1550 -13.541
## Year15                               -0.4132     0.1565 47.2363  -2.640
## Year16                               -0.7093     0.1688 50.6884  -4.202
## D14CommunityMixed                    -0.7146     0.3265 73.0936  -2.189
## D14CommunityDurusdinium-only         -0.6296     0.3928 72.9975  -1.603
## CommunityMixed                        0.7009     0.2845 70.5547   2.464
## CommunityDurusdinium-only             0.7527     0.3508 70.9827   2.146
## Year15:D14CommunityMixed              1.5019     0.3806 64.8657   3.946
## Year16:D14CommunityMixed              0.8259     0.2396 53.9251   3.447
## Year15:D14CommunityDurusdinium-only   1.2674     0.5110 61.4046   2.480
## Year16:D14CommunityDurusdinium-only   0.5426     0.2368 48.9752   2.292
## Year15:CommunityMixed                -0.7379     0.3653 70.6556  -2.020
## Year15:CommunityDurusdinium-only     -0.8090     0.4572 64.7392  -1.770
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## Year15                              0.011206 *  
## Year16                              0.000107 ***
## D14CommunityMixed                   0.031805 *  
## D14CommunityDurusdinium-only        0.113266    
## CommunityMixed                      0.016178 *  
## CommunityDurusdinium-only           0.035332 *  
## Year15:D14CommunityMixed            0.000198 ***
## Year16:D14CommunityMixed            0.001108 ** 
## Year15:D14CommunityDurusdinium-only 0.015876 *  
## Year16:D14CommunityDurusdinium-only 0.026258 *  
## Year15:CommunityMixed               0.047162 *  
## Year15:CommunityDurusdinium-only    0.081500 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 14 columns / coefficients
    #step(SH_Communities)
    
    SH.emmc<-emmeans(SH_Communities, ~Year*D14Community*Community)
    SH_groups<-cld(SH.emmc)
    SH_groups<-as.data.frame(SH_groups[complete.cases(SH_groups),])
    SH_groups<-SH_groups[order(SH_groups$Year,
                               SH_groups$D14Community,
                               SH_groups$Community),]
    SH_groups
    #write.csv(SH_groups,"SH_Communities_multicomp.csv")

SH_Comties <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
    
      facet_grid(~D14Community, labeller = labeller(D14Community=c(
                      `Cladocopium-only` = "Cladocopium-only (2014)",
                      `Mixed` = "Mixed (2014)", 
                      `Durusdinium-only` = "Durusdinium-only (2014)"))) +
        
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.8),
            name="\n Total symbiont to host cell ratio (S/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(14, 15, 16),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) 
SH_Comties

#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)

8a. Changes in Cladocopium abundance

CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(DominantCommunity))) + ggthe_bw+ ggtitle("a.")+
  #geom_point()+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3","purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Cladocopium", "Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.8),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.7),
            name="\n Cladocopium to host cell ratio (C/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" ", " ", " ")) 
#CH_Com

CH_Com2<-CH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
                   `Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
                   `Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))


CH<-lmer(logC.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
    summary(CH)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * D14Dominant + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 113.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.70311 -0.36865 -0.09556  0.50582  1.89250 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.1182   0.3438  
##  Residual             0.4032   0.6350  
## Number of obs: 55, groups:  Colony, 21
## 
## Fixed effects:
##                                         Estimate Std. Error       df t value
## (Intercept)                             -1.54068    0.23810 44.78323  -6.471
## Year15                                  -0.25625    0.28031 26.58185  -0.914
## Year16                                  -0.73614    0.32339 24.71137  -2.276
## CommunityMixed                          -0.32887    0.37022 44.99572  -0.888
## D14DominantDurusdinium-dominated        -1.09559    0.41353 43.96389  -2.649
## Year15:CommunityMixed                    0.59388    0.61658 32.70059   0.963
## Year16:CommunityMixed                   -0.03386    0.50241 28.81910  -0.067
## Year15:D14DominantDurusdinium-dominated -0.54917    0.68606 30.91184  -0.800
## Year16:D14DominantDurusdinium-dominated  0.88798    0.66184 31.43058   1.342
##                                         Pr(>|t|)    
## (Intercept)                             6.36e-08 ***
## Year15                                    0.3688    
## Year16                                    0.0318 *  
## CommunityMixed                            0.3791    
## D14DominantDurusdinium-dominated          0.0112 *  
## Year15:CommunityMixed                     0.3425    
## Year16:CommunityMixed                     0.9467    
## Year15:D14DominantDurusdinium-dominated   0.4296    
## Year16:D14DominantDurusdinium-dominated   0.1893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmntM D14DD- Y15:CM Y16:CM Y15:D1
## Year15      -0.715                                                 
## Year16      -0.579  0.492                                          
## CommuntyMxd -0.622  0.502  0.369                                   
## D14DmnntDr- -0.019 -0.038  0.003 -0.537                            
## Yr15:CmmntM  0.318 -0.474 -0.224 -0.512  0.275                     
## Yr16:CmmntM  0.400 -0.331 -0.659 -0.643  0.345  0.357              
## Yr15:D14DD-  0.006  0.017  0.000  0.255 -0.468 -0.705 -0.185       
## Yr16:D14DD- -0.021  0.011  0.011  0.308 -0.509 -0.162 -0.437  0.306
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
    #step(CH)

    CH.emmc<-emmeans(CH, ~Year*Community* D14Dominant)
    CH_groups<-cld(CH.emmc)
    CH_groups<-as.data.frame(CH_groups[complete.cases(CH_groups),])
    CH_groups<-CH_groups[order(CH_groups$Year,CH_groups$D14Dominant, CH_groups$Community),]
    CH_groups
    #write.csv(CH_groups,"CH_groups.csv")

9a. Changes in Durusdinium abundance

DH_Com <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) + 
      ggthe_bw+ ggtitle("b.")+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Detected symbionts",
                        labels=c("Mixed", "Durusdinium"))+
      theme(legend.position=c(0.82,0.2),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-6, -0.7),
            name="\n Durusdinium to host cell ratio (D/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
                           
# DH_Com

DH_Com2<-DH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
                   `Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
                   `Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))


DH_LM<-lmer(logD.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
    summary(DH_LM)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * D14Dominant + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 66.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0953 -0.5756  0.0079  0.3844  3.6310 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.0000   0.0000  
##  Residual             0.1687   0.4107  
## Number of obs: 58, groups:  Colony, 22
## 
## Fixed effects:
##                                                            Estimate Std. Error
## (Intercept)                                                -4.86573    0.16766
## Year15                                                      0.07049    0.33532
## Year16                                                      3.21459    0.23711
## CommunityDurusdinium-only                                   0.27806    0.50992
## D14DominantDurusdinium-dominated                            3.15614    0.23711
## Year15:CommunityDurusdinium-only                           -0.35506    0.32897
## Year16:CommunityDurusdinium-only                           -0.22094    0.38416
## Year15:D14DominantDurusdinium-dominated                     0.36688    0.42745
## Year16:D14DominantDurusdinium-dominated                    -3.09697    0.41068
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated -0.08557    0.45916
##                                                                  df t value
## (Intercept)                                                48.00000 -29.021
## Year15                                                     48.00000   0.210
## Year16                                                     48.00000  13.557
## CommunityDurusdinium-only                                  48.00000   0.545
## D14DominantDurusdinium-dominated                           48.00000  13.311
## Year15:CommunityDurusdinium-only                           48.00000  -1.079
## Year16:CommunityDurusdinium-only                           48.00000  -0.575
## Year15:D14DominantDurusdinium-dominated                    48.00000   0.858
## Year16:D14DominantDurusdinium-dominated                    48.00000  -7.541
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 48.00000  -0.186
##                                                            Pr(>|t|)    
## (Intercept)                                                 < 2e-16 ***
## Year15                                                        0.834    
## Year16                                                      < 2e-16 ***
## CommunityDurusdinium-only                                     0.588    
## D14DominantDurusdinium-dominated                            < 2e-16 ***
## Year15:CommunityDurusdinium-only                              0.286    
## Year16:CommunityDurusdinium-only                              0.568    
## Year15:D14DominantDurusdinium-dominated                       0.395    
## Year16:D14DominantDurusdinium-dominated                    1.09e-09 ***
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated    0.853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 CmmnD- D14DD- Y15:CD Y16:CD Y15:D1 Y16:D1
## Year15      -0.500                                                        
## Year16      -0.707  0.354                                                 
## CmmntyDrsd-  0.000  0.000 -0.232                                          
## D14DmnntDr- -0.707  0.354  0.500 -0.232                                   
## Yr15:CmmnD-  0.000  0.000  0.000 -0.293  0.360                            
## Yr16:CmmnD-  0.000  0.000  0.000 -0.753  0.309  0.389                     
## Yr15:D14DD-  0.392 -0.784 -0.277  0.129 -0.555 -0.500 -0.171              
## Yr16:D14DD-  0.408 -0.204 -0.577  0.671 -0.577 -0.208 -0.713  0.320       
## CmD-:D14DD-  0.000  0.000  0.258 -0.900  0.000  0.000  0.558  0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
    #step(DH_LM)

    DH.emmc<-emmeans(DH_LM, ~Year*Community* D14Dominant)
    DH_groups<-cld(DH.emmc)
    DH_groups<-as.data.frame(DH_groups[complete.cases(DH_groups),])
    DH_groups<-DH_groups[order(DH_groups$Year,DH_groups$D14Dominant, DH_groups$Community),]
    DH_groups
    #write.csv(DH_groups,"DH_groups.csv")

Figure 4: Changes in C/H and D/H

Figure_4<-grid.arrange(CH_Com2, DH_Com2,nrow=2)

# ggsave(file="Outputs/Figure_4.svg", plot=Figure_4, width=7, height=4)

Symbiodiiaceae ITS2 type + ORF lineage data analysis

Added Pocillopora ORF lineage and Symbiodiniaceae ITS2 type hosted in 2014

  #data.Uva$Counts<-paste(data.Uva$DominantCommunity, data.Uva$ORFCommunity, sep = "_")
  data.Uva$ORFCommunity<-paste(data.Uva$ORF, data.Uva$Community, sep = "_")
  
  data.Uva$ORFCommunity<-factor(data.Uva$ORFCommunity,
                                levels = c("Type_3_Cladocopium-only", 
                                           "Type_1_Cladocopium-only",
                                           "Type_1_Mixed", 
                                           "Type-1_Durusdinium-only"))
  
  data.Uva$ORFCommunity2<-paste(data.Uva$ORF,
                              data.Uva$DominantCommunity, sep = "_")
  data.Uva$ORFCommunity2<-as.factor(data.Uva$ORFCommunity2)
 
  data.Uva$DominantITS2[which(
    (data.Uva$DominantCommunity=="Durusdinium-dominated") &
    (data.Uva$ORF=="Type_1")
      )] <- "D1-dominated"
  
  data.Uva$DominantITS2[which(
    (data.Uva$DominantCommunity=="Cladocopium-dominated") &
    (data.Uva$ORF=="Type_1")
      )] <- "C1bc-dominated"
  
    data.Uva$DominantITS2[which(
    (data.Uva$DominantCommunity=="Cladocopium-dominated") &
    (data.Uva$ORF=="Type_3")
      )] <- "C1d-dominated"
  
   data.Uva$DominantITS2<-as.factor(data.Uva$DominantITS2)
   
  Summary_Community4<-data.Uva %>%
        group_by(Year.2, DominantITS2) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community4 
  Summary_Community5<-data.Uva %>%
        group_by(ORF, Year.2, DominantITS2) %>%
        dplyr::summarise(total.count=n()) 
  Summary_Community5

5b. Changes in the number of colonies dominated by each genus

  • Fisher test including initially (Aug 2014) 50/50 colonies. Both colonies placed in the D-dominated category to be conservative.
# Only type 1 colonies (n=23), p-value = 0.03514

Input =("
Year    C     D
2014    9    14
2016    2    21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.03514
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.105302 70.563784
## sample estimates:
## odds ratio 
##     6.4764
# Type 1 and type 3 colonies: p-value = 0.1064

Input =("
Year    C     D
2014    15    14
2016    8     21
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.8336155 9.7658946
## sample estimates:
## odds ratio 
##   2.760908

6b. Changes in colony proportions based on ITS2 profiles (2014) and genus dominance

  • Histograms based on ITS2 detected symbionts in 2014
data.Uva$ITS2_14<-factor(data.Uva$ITS2_14, levels = c("C1d", "C1b-c", "D1"))
data.Uva$ITS2_16<-factor(data.Uva$ITS2_16, levels = c("C1d", "C1b-c", "D1"))

Histo_Dprop_ITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(ITS2_14)), shape=(D14Community))  + ggthe_bw +
      geom_histogram(binwidth = 0.05, boundary=0) +
      facet_grid(~Year.2, labeller = labeller(Year.2=c(
                    `2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
      scale_fill_manual(values=c("blue3", "purple", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "Initial ITS2 community (August 2014)", 
                        labels=c("Type 3 + C1d", "Type 1 + C1b-c", "Type 1 + D1")) +
      scale_x_continuous(name= expression(
                        italic(Durusdinium)~proportion~(D/H)/(S/H)),
                        breaks=(seq(0,1,by=0.1))) +
      scale_y_continuous((name= "Number of colonies"),
                         breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
Histo_Dprop_ITS2 +coord_flip()

  • Histograms based on ORF lineage
Histo_Dprop_ORF <- ggplot(data.Uva, aes (D.Prp , fill=(ORF)), shape=(D14Community))  + ggthe_bw +
      geom_histogram(binwidth = 0.05, boundary=0) +
      facet_grid(~Year.2, labeller = labeller(Year.2=c(
                    `2014` = "August 2014", 
                    `2015` = "August 2015", `2016` = "April 2016"))) +
      scale_fill_manual(values=c("black","gray"),
                        guide = guide_legend(title.position = "top", nrow = 1),
                        name = "ORF type", 
                        labels=c("Type 3", "Type 1")) +
      scale_x_continuous(name= expression(
                        italic(Durusdinium)~proportion~(D/H)/(S/H)),
                        breaks=(seq(0,1,by=0.1))) +
      scale_y_continuous((name= "Number of colonies"),
                         breaks=seq(0,15,by=5)) +
      theme(legend.position="bottom",
            legend.title.align=0.5,
            strip.background =element_rect(fill=NA))
Histo_Dprop_ORF +coord_flip()

  • Durusdinium proportion trajectory
# By ITS2 type
Traject_Proportion_FacetITS2 <- ggplot (data.Uva) +
   ggthe_bw + 
    geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+ 
    facet_wrap(~ITS2_14) +
    geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3, 
              height = 0.05, width = 0.15) +
    #labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
    theme(legend.position="none")
Traject_Proportion_FacetITS2

# By Depth
Traject_Proportion_FacetORF <- ggplot (data.Uva) +    ggthe_bw +
      geom_line(aes (x=Year.2, D.Prp, group=Colony, linetype=(ORF)))+
      geom_jitter(aes(Year.2, D.Prp, colour=(Depht), shape=ORF),
              height = 0.00, width = 0.25, size=3, alpha=0.7) +  
      scale_colour_gradient(low = "red", high = "blue")
Traject_Proportion_FacetORF

Figure XX: ITS2

Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ITS2_14))) +
        ggthe_bw + ggtitle("a.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Dominant ITS2 community (August 2014)", 
                          labels=c("C1d", "C1b-c", "D.glynnii")) +
        scale_x_continuous(name= expression(
                          italic(Durusdinium)~proportion),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2 

Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ITS2_16))) +
        ggthe_bw + ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Dominant ITS2 community (April 2016)", 
                          labels=c("C1d", "C1b-c", "D.glynnii")) +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2

Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , 
                                fill=ITS2_14, shape=Mortality, group=(Colony)))+
        ggthe_bw + ggtitle("b.") +
        geom_line(aes(colour=ITS2_14), linetype=2) +
        geom_jitter(height = 0, width = 0.25, alpha=0.4) +
        
        scale_colour_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Initial ITS2 community (August 2014)", 
                          labels=c("C1d", "C1b-c", "D.glynnii")) +
        
        scale_fill_manual(values=c("blue3","purple", "red3"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Initial ITS2 community (August 2014)", 
                          labels=c("C1d", "C1b-c", "D.glynnii")) +
        scale_shape_manual(values=c(21, 22, 24),
                      labels=c("No", "Partial", "Almost total"))+
        
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.01,0.01)) +
        
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Traject_ProportionITS2


Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(2/9, 5/9, 2/9))

# ggsave(file="Outputs/Figure_ITS2.svg", plot=Figure_2_ITS2, width=6.5, height=5)
Histo2014_DpropITS2b<- Histo2014_DpropITS2 + facet_grid(ORF~.)

Histo2016_DpropITS2b<- Histo2016_DpropITS2 + facet_grid(ORF~.)

Traject_ProportionITS2b<- Traject_ProportionITS2 + facet_grid(ORF~.)

Figure_2_ITS2b<-grid.arrange(Histo2014_DpropITS2b, Traject_ProportionITS2b,Histo2016_DpropITS2b,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_2_ITS2b.svg", plot=Figure_2_ITS2b, width=7, height=4)

Figure XX: ORF

Histo2014_DpropORF <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ORFCommunity))) +
        ggthe_bw + coord_flip()+ ggtitle("a.")+
        geom_histogram(binwidth = 0.05, boundary=0, colour="black") +
        scale_fill_manual(values=c("gray", "white","snow", "mintcream"),
                       guide = guide_legend(title.position = "top", nrow = 2),
                      name = "")+ 
                       #labels=c("Type 3", "Type 1")) +
        scale_x_continuous(#name= expression(
                          #italic(Durusdinium)~proportion~((D/H)/(S/H))),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                          limits = c(0,16),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA),
              panel.border = element_blank(),
              axis.line = element_line(colour = "black"))
#Histo2014_DpropORF 

Histo2016_DpropORF <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ORFCommunity))) +
        ggthe_bw + coord_flip()+ ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0, 
                       colour="black") +
        scale_fill_manual(values=c("gray", "white","snow", "mintcream"),
                       guide = guide_legend(title.position = "top", nrow = 2),
                      name = "")+ 
                       #labels=c("Type 3", "Type 1")) +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,16),
                          breaks=seq(0,15,by=5)) +
          theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA),
              panel.border = element_blank(),
              axis.line = element_line(colour = "black"))
#Histo2016_DpropORF

Traject_ProportionORF_a<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp))+
        ggthe_bw + ggtitle("b.") +
        geom_line(aes(colour=(ORF), group=(Colony)), linetype=2) +
        scale_colour_manual(values=c("black", "gray"),
                          guide = guide_legend(title.position = "top", nrow = 2),
                          name = "Pocillopora type", 
                          labels=c("Type 3", "Type 1")) +
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.01,0.01)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
  
  Traject_ProportionORF <- Traject_ProportionORF_a + 
  geom_jitter((aes(colour=factor(ORF))), 
                    height = 0, width = 0.3, color="black", size=1, shape=21) +
        scale_fill_manual(values=c("white", "gray"))
  
   #Traject_ProportionORF

Figure_2_ORF<-grid.arrange(Histo2014_DpropORF, Traject_ProportionORF, Histo2016_DpropORF,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_2_ORF.svg",plot=Figure_2_ORF, width=7, height=3.8, dpi=300)

Figure XX: ITS2 + Moratality

data.Uva$Mor_ITS2<-paste(data.Uva$DominantITS2, data.Uva$Mortality, sep = "_")

data.Uva$Mor_ITS2<-factor(data.Uva$Mor_ITS2, levels = c(
  "C1d-dominated_No", "C1d-dominated_Partial", "C1d-dominated_Almost total",  
  "C1bc-dominated_No", "C1bc-dominated_Partial",  
  "D1-dominated_No", "D1-dominated_Partial", "D1-dominated_Almost total"))

summary(data.Uva$Mor_ITS2)
##           C1d-dominated_No      C1d-dominated_Partial 
##                         11                          3 
## C1d-dominated_Almost total          C1bc-dominated_No 
##                          4                         17 
##     C1bc-dominated_Partial            D1-dominated_No 
##                          3                         43 
##       D1-dominated_Partial  D1-dominated_Almost total 
##                          3                          3
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + coord_flip()+ ggtitle("a.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("skyblue1", 
                          "darkorchid1",
                          "red1"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (August 2014)") +
      
        scale_x_continuous(name= expression(
                          italic(Durusdinium)~proportion~((D/H)/(S/H))),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2 

Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + coord_flip()+ ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
       scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
                          "darkorchid3",
                          "red1", "red3", "red4"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2


HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + coord_flip()+ ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
       scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
                          "darkorchid1","darkorchid3",
                          "red1", "red3", "red4"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#HistoAll_DpropITS2


Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
        ggthe_bw + ggtitle("b.") +
        geom_line(aes(colour=ITS2_14),linetype=2) +
        geom_jitter(aes(colour=Mor_ITS2, shape=Mortality),
                    height = 0.00, width = 0.25, alpha=0.9) +
        
        scale_colour_manual(values=c(
                          "darkorchid1", "darkorchid1","darkorchid3",
                          "skyblue1","blue4", "skyblue1", "blue2",
                          "red1", "red4", "red1", "red3" ),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
  
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.01,0.01)) +
        
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))
#Traject_ProportionITS2


Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(3/9, 3/9, 3/9))

# ggsave(file="Outputs/Figure_ITS2_Mortality.svg", plot=Figure_2_ITS2, width=7, height=4)

Figure 3b and 3c: ITS2 + Moratality

data.Uva$Mor_ITS2<-paste(data.Uva$DominantITS2, data.Uva$Mortality, sep = "_")

data.Uva$Mor_ITS2<-factor(data.Uva$Mor_ITS2, levels = c(
   "C1d-dominated_Almost total", "C1d-dominated_Partial","C1d-dominated_No",  
  "C1bc-dominated_Partial","C1bc-dominated_No",
  "D1-dominated_Almost total", "D1-dominated_Partial", "D1-dominated_No" ))

summary(data.Uva$Mor_ITS2)
## C1d-dominated_Almost total      C1d-dominated_Partial 
##                          4                          3 
##           C1d-dominated_No     C1bc-dominated_Partial 
##                         11                          3 
##          C1bc-dominated_No  D1-dominated_Almost total 
##                         17                          3 
##       D1-dominated_Partial            D1-dominated_No 
##                          3                         43
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + ggtitle("a.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
        scale_fill_manual(values=c("black", 
                          "black",
                          "black"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (August 2014)") +
      
        scale_x_continuous(name= expression(
                          italic(Durusdinium)~proportion~((D/H)/(S/H))),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="none",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) + facet_grid(ORF~.)
#Histo2014_DpropITS2 

Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
       scale_fill_manual(values=c("white", "gray", "black",
                          "gray",
                          "white", "gray", "black"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="none",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) +  facet_grid(ORF~.)
#Histo2016_DpropITS2


HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0) +
       scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
                          "darkorchid1","darkorchid3",
                          "red1", "red3", "red4"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#HistoAll_DpropITS2

HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
        ggthe_bw + ggtitle("c.")+
        geom_histogram(binwidth = 0.05, boundary=0, color = "black") +
       scale_fill_manual(values=c("white", "gray", "black",
                          "gray","black",
                          "white", "gray", "black"),
                          guide = guide_legend(title.position = "top", nrow = 3),
                          name = "Dominant ITS2 community (April 2016)") +
        scale_x_continuous(name="",
                           breaks=(seq(0,1,by=0.1)),
                           expand=c(0.01,0.01)) +    
        scale_y_continuous((name= "Number of colonies"),
                          limits = c(0,14),
                          breaks=seq(0,14,by=2)) +
        theme(legend.position="bottom",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) + 
        facet_grid(ORF~Year) + theme(legend.position = "none")
#HistoAll_DpropITS2
#ggsave(file="Outputs/Panels.svg", plot=HistoAll_DpropITS2, width=6, height=2.8)



Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
        ggthe_bw + ggtitle("b.") +
        geom_line(linetype=2) +
        geom_jitter(aes(colour=Mortality, shape=ORF), size=1.5,
                    height = 0.00, width = 0.25, alpha=1) +
        
        scale_colour_manual(values=c("black", "gray","white")) +
  
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.02,0.02)) +
        
        theme(legend.position="none",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))# + facet_grid(ORF~.)
#Traject_ProportionITS2
#ggsave(file="Outputs/Trajectory.svg", plot=Traject_ProportionITS2, width=6.5, height=3)

Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(2.5/9, 4/9, 2.5/9))

# ggsave(file="Outputs/Figure_ITS2_Mortality2.svg", plot=Figure_2_ITS2, width=7, height=4)
DpropORF <- ggplot(data.Uva, aes (D.Prp , fill=(Mortality))) +
        ggthe_bw + ggtitle("a.")+
        geom_histogram(binwidth = 0.1, boundary=0, color="black") +
        scale_fill_manual(values=c("black", 
                          "gray",
                          "white"),
                          guide = guide_legend(title.position = "top", nrow = 1),
                          name = "Colony condition") +
      
        scale_x_continuous(name= expression(
                          italic(Durusdinium)~proportion~((D/H)/(S/H))),
                          breaks=(seq(0,1,by=0.1)),
                          expand=c(0.01,0.01)) +
        scale_y_continuous((name= "Number of colonies 
                            \n (Aug 2014) 
                            \n (Aug 2014) \n (Aug 2014)"),
                          limits = c(0,17),
                          breaks=seq(0,15,by=5)) +
        theme(legend.position="top",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#DpropORF 

Traject_ORF<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
        ggthe_bw + ggtitle("b.") +
        geom_line(linetype=2) +
        geom_jitter(aes(colour=Mortality, shape=ORF), size=2,
                    height = 0.00, width = 0.25, alpha=1) +
        
        scale_colour_manual(values=c("black", "gray","white")) +
  
        scale_x_continuous(name="", 
                           breaks = c(2014, 2015, 2016),
                           labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
        scale_y_continuous((name= ""),
                           breaks=seq(0,1,by=0.1),
                           expand=c(0.02,0.02)) +
        
        theme(legend.position="none",
              legend.title.align=0.5,
              strip.background =element_rect(fill=NA))# + facet_grid(ORF~.)
#Traject_ORF

Figure_3bc_ORF<-grid.arrange(Traject_ORF, DpropORF,
nrow=2, heights=c(1.2/3, 1.8/3))

#ggsave(file="Outputs/Fig3_v8.svg", plot=Figure_3bc_ORF, width=6.5, height=7)
# Summary dominant genera
  Community_prescence_absence<-data.Uva %>%
        group_by(ORF, Community, Year.2) %>%
        dplyr::summarise(total.count=n()) 
  Community_prescence_absence
  Community_dominance<-data.Uva %>%
        group_by(ORF, DominantCommunity, Year.2) %>%
        dplyr::summarise(total.count=n()) 
  Community_dominance
  Community_dominance<-data.Uva %>%
        group_by(ORF, Year.2, DominantCommunity,
                Community) %>%
        dplyr::summarise(total.count=n()) 
  Community_dominance
  # SH_mean<-data.Uva %>%
  #       group_by(ORF, Year.2, DominantCommunity) %>%
  #       dplyr::summarise(mean=mean(tot.SH)) 
  # SH_mean
  # 
  # SH_SE<-data.Uva %>%
  #       group_by(ORF, Year.2, DominantCommunity) %>%
  #       dplyr::summarise(sd=sd(tot.SH)) 
  # SH_SE

7b. Changes in total Symbiodiniaceae abundance (bleaching): qPCR + ORF

Figure 3: ORF + Dominant symbiont

#SH_2014CoI<-lm(logTot.SH ~ Year * ORF * DominantCommunity,             data=data.Uva)
SH_2014CoI<-lmerTest::lmer(logTot.SH ~ Year * ORF * DominantCommunity + (1|Colony), data=data.Uva)

    summary(SH_2014CoI)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ORF * DominantCommunity + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 75.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3367 -0.5185 -0.0403  0.3713  3.4413 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.01992  0.1411  
##  Residual             0.10475  0.3236  
## Number of obs: 87, groups:  Colony, 29
## 
## Fixed effects:
##                                               Estimate Std. Error       df
## (Intercept)                                   -1.68521    0.11691 77.39879
## Year15                                         0.22667    0.15257 51.72506
## Year16                                        -0.05117    0.26368 68.09623
## ORFType_3                                      0.02507    0.18560 75.84196
## DominantCommunityDurusdinium-dominated         0.14102    0.14919 77.99999
## Year15:ORFType_3                              -0.79408    0.24123 51.72506
## Year16:ORFType_3                              -0.72686    0.32318 62.79186
## Year15:DominantCommunityDurusdinium-dominated -0.06059    0.19555 51.72506
## Year16:DominantCommunityDurusdinium-dominated  0.01196    0.28860 68.27239
##                                               t value Pr(>|t|)    
## (Intercept)                                   -14.415   <2e-16 ***
## Year15                                          1.486   0.1434    
## Year16                                         -0.194   0.8467    
## ORFType_3                                       0.135   0.8929    
## DominantCommunityDurusdinium-dominated          0.945   0.3475    
## Year15:ORFType_3                               -3.292   0.0018 ** 
## Year16:ORFType_3                               -2.249   0.0280 *  
## Year15:DominantCommunityDurusdinium-dominated  -0.310   0.7579    
## Year16:DominantCommunityDurusdinium-dominated   0.041   0.9671    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 ORFT_3 DmnCD- Y15:OR Y16:OR Y15:DC
## Year15      -0.653                                                 
## Year16      -0.373  0.289                                          
## ORFType_3   -0.630  0.411  0.235                                   
## DmnntCmmnD- -0.777  0.511  0.290  0.489                            
## Yr15:ORFT_3  0.413 -0.632 -0.183 -0.650 -0.323                     
## Yr16:ORFT_3  0.305 -0.236 -0.816 -0.483 -0.237  0.373              
## Yr15:DmnCD-  0.509 -0.780 -0.226 -0.321 -0.655  0.493  0.184       
## Yr16:DmnCD-  0.359 -0.264 -0.920 -0.226 -0.463  0.167  0.750  0.339
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
    anova(SH_2014CoI)
    ranova(SH_2014CoI)
    step(SH_2014CoI, reduce.random=FALSE)
## Backward reduced random-effect table:
## 
##              Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                    11 -37.966 97.933                     
## (1 | Colony)          0   10 -38.884 97.768 1.8355  1     0.1755
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                            Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Year:ORF:DominantCommunity          0                                          
## 
## Model found:
## logTot.SH ~ Year * ORF * DominantCommunity + (1 | Colony)
    SH_2014CoI<-lm(logTot.SH ~ Year * ORF * DominantCommunity, data=data.Uva)
    
    SH_I.emmc<-emmeans(SH_2014CoI, ~Year*ORF*DominantCommunity)
    SH_I_groups<-multcomp::cld(SH_I.emmc)
    SH_I_groups<-as.data.frame(SH_I_groups[complete.cases(SH_I_groups),])
    SH_I_groups<-SH_I_groups[order(
      SH_I_groups$Year, SH_I_groups$ORF),]
    SH_I_groups
    #write.csv(SH_I_groups,"Outputs/SH_ORF_Dominant_groups2.csv")
    
# plot(emmeans(SH_2014CoI, ~Year|~DominantCommunity), comparisons = TRUE) + coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~DominantITS2)+ ggthe_bw
SH_2014ComI <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(DominantITS2))) +  ggthe_bw + facet_wrap(~ORF) + ggtitle("a.")+
  #geom_point(alpha=0.5)+
  #geom_line(aes(group=(Colony)))+
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2,
                   position=position_dodge(width=0.5))+
  stat_summary(fun.data = "mean_cl_boot",geom = "point", 
                   position=position_dodge(width=0.5))+
  stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.5)) +
   scale_colour_manual(values=c("blue3","blue3" , "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Dominant symbiont",
                        labels=c("Cladocopium", "Durusdinium"))+
      theme(legend.position="mone",
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-3, -0.1),
            name="\n Total symbiont to host cell ratio (S/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" ", "", " ")) 
SH_2014ComI

#ggsave(file="Figure_3.svg", plot=SH_ITS2_Comunities, width=6, height=4, dpi=300)

8b. Changes in Cladocopium abundance

CH_Com2ITS2 <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(D14Dominant))) + ggthe_bw+ ggtitle("c.") + facet_grid(~ORF) + #geom_point(alpha=0.5)+
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue3", "red3"),
                         guide = guide_legend(title.position = "top", nrow = 3),
                         name = "Dominant symbiont",
                         labels=c("Cladocopium", "Durusdinium"))+
     theme(legend.position="mone",
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-6, -0.1),
            name="\n Cladocopium to host cell ratio (C/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                         labels = c(" ", " ", " ")) 
CH_Com2ITS2

CH_I<-lmerTest::lmer(logC.SH ~ Year* DominantITS2 + (1|Colony), data=data.Uva)
    summary(CH_I)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * DominantITS2 + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 108.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.89915 -0.40700 -0.07625  0.38592  1.92091 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.06891  0.2625  
##  Residual             0.38254  0.6185  
## Number of obs: 55, groups:  Colony, 21
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                      -1.68510    0.22354 44.93743  -7.538 1.65e-09
## Year15                            0.22932    0.29156 28.78344   0.787 0.437988
## Year16                           -0.12883    0.50336 39.26788  -0.256 0.799334
## DominantITS2C1d-dominated         0.02496    0.35386 44.57426   0.071 0.944077
## DominantITS2D1-dominated         -1.27282    0.35244 45.62565  -3.611 0.000755
## Year15:DominantITS2C1d-dominated -0.79673    0.46100 28.78344  -1.728 0.094658
## Year16:DominantITS2C1d-dominated -0.64921    0.61716 35.78277  -1.052 0.299886
## Year15:DominantITS2D1-dominated  -0.44588    0.49919 31.70568  -0.893 0.378478
## Year16:DominantITS2D1-dominated   0.31093    0.62644 42.90497   0.496 0.622184
##                                     
## (Intercept)                      ***
## Year15                              
## Year16                              
## DominantITS2C1d-dominated           
## DominantITS2D1-dominated         ***
## Year15:DominantITS2C1d-dominated .  
## Year16:DominantITS2C1d-dominated    
## Year15:DominantITS2D1-dominated     
## Year16:DominantITS2D1-dominated     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 DITS2C DITS2D Y15:DITS2C Y16:DITS2C Y15:DITS2D
## Year15      -0.652                                                             
## Year16      -0.376  0.290                                                      
## DmnnITS2C1- -0.632  0.412  0.238                                               
## DmnnITS2D1- -0.631  0.414  0.238  0.398                                        
## Y15:DITS2C1  0.412 -0.632 -0.183 -0.651 -0.262                                 
## Y16:DITS2C1  0.307 -0.236 -0.816 -0.486 -0.194  0.373                          
## Y15:DITS2D1  0.382 -0.584 -0.170 -0.241 -0.606  0.369      0.138               
## Y16:DITS2D1  0.339 -0.233 -0.816 -0.214 -0.538  0.147      0.665      0.345
    anova (CH_I)
    #step(CH_I)

    CH_I.emmc<-emmeans(CH_I, ~Year*DominantITS2)
    CH_I_groups<-cld(CH_I.emmc)
    CH_I_groups<-as.data.frame(CH_I_groups[complete.cases(CH_I_groups),])
    CH_I_groups<-CH_I_groups[order(CH_I_groups$DominantITS2,CH_I_groups$Year),]
    CH_I_groups
    #write.csv(CH_I_groups,"CH_I_groups.csv")

9b. Changes in Durusdinium abundance

DH_Com2ITS2 <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(D14Dominant))) +  # geom_point(alpha=0.5)+
      ggthe_bw+ ggtitle("d.")+ facet_grid(~ORF) +
    
      stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
      stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3))+
      stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
      
      scale_colour_manual(values=c("blue", "red3"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Dominant symbionts",
                        labels=c("Cladocopium-dominated", "Durusdinium-dominated"))+
      theme(legend.position=c(0.75,0.5),
            legend.box.background = element_rect(),
            strip.background =element_rect(fill="white"))+
   
      scale_y_continuous(limits = c(-6, -0.1),
            name="\n Durusdinium to host cell ratio (D/H)", 
            labels = (math_format(10^.x)),
            expand = c(0.0, 0)) + 
      scale_x_continuous(name="", 
                         breaks = c(2014, 2015, 2016),
                        # labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
                        labels = c(" ", " ", " "))

DH_Com2ITS2

DH_LM_ITS2<-lmer(logD.SH ~ Year* DominantITS2 + (1|Colony), data=data.Uva)
    summary(DH_LM_ITS2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * DominantITS2 + (1 | Colony)
##    Data: data.Uva
## 
## REML criterion at convergence: 62.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1646 -0.5042  0.0345  0.3558  3.7511 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Colony   (Intercept) 0.000    0.0000  
##  Residual             0.158    0.3975  
## Number of obs: 58, groups:  Colony, 22
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                     -4.86573    0.16229 52.00000 -29.981  < 2e-16
## Year15                           0.07049    0.32459 52.00000   0.217    0.829
## Year16                           2.95052    0.42939 52.00000   6.871 7.93e-09
## DominantITS2D1-dominated         3.26614    0.19398 52.00000  16.838  < 2e-16
## Year15:DominantITS2D1-dominated  0.14077    0.35768 52.00000   0.394    0.696
## Year16:DominantITS2D1-dominated -2.96086    0.45076 52.00000  -6.569 2.41e-08
##                                    
## (Intercept)                     ***
## Year15                             
## Year16                          ***
## DominantITS2D1-dominated        ***
## Year15:DominantITS2D1-dominated    
## Year16:DominantITS2D1-dominated ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year15 Year16 DITS2D Y15:DI
## Year15      -0.500                            
## Year16      -0.378  0.189                     
## DmnnITS2D1- -0.837  0.418  0.316              
## Y15:DITS2D1  0.454 -0.907 -0.171 -0.542       
## Y16:DITS2D1  0.360 -0.180 -0.953 -0.430  0.233
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
    #step(DH_LM_ITS2)

    DH_I.emmc<-emmeans(DH_LM_ITS2, ~Year*DominantITS2)
    DH_groups_I<-cld(DH_I.emmc)
    DH_groups_I<-as.data.frame(DH_groups_I[complete.cases(DH_groups_I),])
    DH_groups_I<-DH_groups_I[order(DH_groups_I$DominantITS2,DH_groups_I$Year),]
    DH_groups_I
    #write.csv(DH_groups_I,"Output/DH_I_groups.csv")

Figure 4: qPCR by Pocillopora lineage

Figure_4<-grid.arrange(SH_2014ComI, CH_Com2ITS2, DH_Com2ITS2,nrow=3)

Figure_4_InitailCommunity<-grid.arrange(CH_Com2ITS2, DH_Com2ITS2,nrow=2)

# ggsave(file="Outputs/Figure_3c.svg", plot=Figure_3, width=7, height=11)

# ggsave(file="Outputs/Figure_3d.png", plot=Figure_3_InitailCommunity, width=7, height=7)

10. Supplmenmtary C/H D/H dynamics

Select data

log_Log.data<-data.Uva[,c("Colony","Year.2", "C.SH", "D.SH", 
                          "D14Community", "Community","ITS2_14",
                          "ORF", "DominantITS2", "Mortality", "Mor_ITS2")]

log_Log.data<-log_Log.data %>%
  filter(Year.2 != "2015")

log_Log.data$Yr<-as.factor(log_Log.data$Year.2)
log_Log.data$Yr_M<-as.factor(
                            paste(log_Log.data$Yr,
                                  log_Log.data$Mortality,
                                  sep = "_"))

log_Log.data$logC<-log10(log_Log.data$C.SH)
log_Log.data$logD<-log10(log_Log.data$D.SH)
log_Log.data$logC[which(log_Log.data$C.SH==0)] <- -6
log_Log.data$logD[which(log_Log.data$D.SH==0)] <- -6
  • ITS2-ORF
Community_Changes1 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14, shape=ORF)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  scale_colour_manual(values=c("blue","purple", "red"),
                        guide = guide_legend(title.position = "top", nrow = 3),
                        name = "Dominant ITS2 types (2014)",
                        labels=c("C1d", "C1b-c", "D1"))+
  
  ggthe_bw + theme(legend.position="bottom", 
                    legend.box = "horizontal")+
  scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
  scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
  
  geom_point(alpha=0.5, size=4) +
  geom_path(size = 1,arrow =
              arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
  
 annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1D : 10,000C", "1D : 1,000C", "1D : 100C",  "1D : 10C"), angle = 45)
  
Community_Changes1

#ggsave(file="FigAndrewsFavoriteGraph_ITS2.png", plot=Community_Changes1, width=5.5, height=7.1)
  • ITS2-Mortality
Community_Changes2 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  scale_colour_manual(values=c("blue","purple", "red"),
        guide = guide_legend(title.position = "top", nrow = 3),
                         name = "Initial ITS2",
                        labels=c("Cd1", "C1b-c", "D1"))+
  #scale_colour_manual(values=c("skyblue1", "blue2", "blue4",
  #                        "darkorchid1","darkorchid3",
  #                        "red1", "red3", "red4"))+
  
  ggthe_bw + 
  theme(legend.position="bottom", 
                   legend.box = "vertical")+
  #theme(legend.position="none")+
  
  scale_x_continuous(limits = c(-6,-0.1), 
                     "Cladocopium abundance (log10 C/H))")+
  scale_y_continuous(limits = c(-6,-0.1), 
                     "Durusdinium abundance (log10 D/H))") +
  geom_point(aes(shape=Mortality),
             alpha=0.5, size=3) +
  geom_path(size = 0.5 , arrow =
             arrow(type = "open", angle = 30,
                  length = unit(0.1, "inches")))+
  
 annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1D : 10,000C", "1D : 1,000C", "1D : 100C",  "1D : 10C"), angle = 45)
  
Community_Changes2

#ggsave(file="LogLog_labels.svg", plot=Community_Changes2, width=5.5, height=6)
  • Colony
Community_Changes3 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=Colony)) +

  geom_abline(intercept=0, slope=1, linetype=1, color="black")+
  geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=1, slope=1, linetype=2, color="black")+
  geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
  geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
  geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
  #scale_colour_manual(values=c("blue","purple", "red"),
   #                     guide = guide_legend(title.position = "top", nrow = 3),
    #                    name = "Detected symbionts (2014)",
     #                   labels=c("Cladocopium", "Mixed", "Durusdinium"))+
  
  ggthe_bw + theme(legend.position="none")+
  scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
  scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
  
  geom_point(aes(shape=Yr), alpha=0.5, size=4) +
  geom_path(size = 1,arrow =
              arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
  
  annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
    "10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
  annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
    "1D : 10,000C", "1D : 1,000C", "1D : 100C",  "1D : 10C"), angle = 45)
  
Community_Changes3

# ggsave(file="LogLog_Colony.svg", plot=Community_Changes3, width=6, height=6)

5c. Changes in the number of colonies dominated by each genus considering mortality

  • Fisher test including initially (Aug 2014) 50/50 colonies. Both colonies placed in the D-dominated category to be conservative.
# Only type 1 colonies (n=23)

Input =("
Year    C     D
2014    9    14
2016    2    18
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.03936
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.9354645 61.0268633
## sample estimates:
## odds ratio 
##   5.560207
# Type 1 and type 3 colonies

Input =("
Year    C     D
2014    15    14
2016    4     18
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.01988
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.145585 23.762126
## sample estimates:
## odds ratio 
##   4.670456
# Spp changes

Input =("
Year    I     II
2014    23    6
2016    20    2
")

Matriz <- as.matrix(read.table(textConnection(Input),
                               header=TRUE, 
                               row.names=1))
fisher.test(Matriz,
            alternative="two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Matriz
## p-value = 0.4399
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.03477518 2.50826671
## sample estimates:
## odds ratio 
##  0.3901317